*******************************************
*****Code for county-level analysis********
*******************************************


use ...county_data.dta, clear 


** Table 1. Estimated impact of the GECP on grassland quality: county and household levels
*** Col. (1) DID 
xi: areg lnNDVI_mn ///
	P_x_T ///
	rain_m5 rain_m6 rain_m7 rain_m8 rain_m9 rain_m10 Tavg_m5 Tavg_m6 Tavg_m7 Tavg_m8 Tavg_m9 Tavg_m10 PSDI_m5 PSDI_m6 PSDI_m7 PSDI_m8 PSDI_m9 PSDI_m10  fiscal_ave ///
	i.year ///
	if (year<=2013 & year>=2008) & (treatgroup2==1|treat==0),a(county_num) rob cluster(county_num) 

*** Col. (2) RE 
xtset county_num
xi: xtreg lnNDVI_mn ///
          P_x_T ///
	      rain_m5 rain_m6 rain_m7 rain_m8 rain_m9 rain_m10 Tavg_m5 Tavg_m6 Tavg_m7 Tavg_m8 Tavg_m9 Tavg_m10 PSDI_m5 PSDI_m6 PSDI_m7 PSDI_m8 PSDI_m9 PSDI_m10 fiscal_ave ///
 	      i.year if (year<=2013 & year>=2008) & (treatgroup2==1|treat==0), re rob cluster(county_num) 


		  
** Table 4. Estimated Impacts of GECP-I on Grazing Scale-County level 
*** col.(1): log(cattle)	
xi: areg lncattle ///
     P_x_T ///
	 rain_m5 rain_m6 rain_m7 rain_m8 rain_m9 rain_m10 Tavg_m5 Tavg_m6 Tavg_m7 Tavg_m8 Tavg_m9 Tavg_m10 PSDI_m5 PSDI_m6 PSDI_m7 PSDI_m8 PSDI_m9 PSDI_m10 fiscal_ave ///
	 i.year ///
	 if (year<=2013 & year>=2008) & (treatgroup2==1|treat==0),a(county_num) rob cluster(county_num) 
	 
*** col.(2): log(sheep)	 
xi: areg lnsheep ///
     P_x_T ///
	 rain_m5 rain_m6 rain_m7 rain_m8 rain_m9 rain_m10 Tavg_m5 Tavg_m6 Tavg_m7 Tavg_m8 Tavg_m9 Tavg_m10 PSDI_m5 PSDI_m6 PSDI_m7 PSDI_m8 PSDI_m9 PSDI_m10 fiscal_ave ///
	 i.year ///
	 if (year<=2013 & year>=2008) & (treatgroup2==1|treat==0),a(county_num) rob cluster(county_num) 

** Table 6. Estimated heterogeneous impacts of GECP on grassland quality: county level
***Panel A. by rural road intensity 
	** Col.1 Low 
   xi: areg lnNDVI_mn ///
    P_x_T ///
	rain_m5 rain_m6 rain_m7 rain_m8 rain_m9 rain_m10 Tavg_m5 Tavg_m6 Tavg_m7 Tavg_m8 Tavg_m9 Tavg_m10 PSDI_m5 PSDI_m6 PSDI_m7 PSDI_m8 PSDI_m9 PSDI_m10 fiscal_ave ///
	i.year   ///
	if (year<=2013 & year>=2008) & (treatgroup2==1|treat==0) & roadv_group==1  & lnNDVI_mn!=., a(county_num) rob cluster(county_num)
	
    ** Col.(2) Medium 
   xi: areg lnNDVI_mn ///
    P_x_T ///
	rain_m5 rain_m6 rain_m7 rain_m8 rain_m9 rain_m10 Tavg_m5 Tavg_m6 Tavg_m7 Tavg_m8 Tavg_m9 Tavg_m10 PSDI_m5 PSDI_m6 PSDI_m7 PSDI_m8 PSDI_m9 PSDI_m10 fiscal_ave ///
	i.year   ///
	if (year<=2013 & year>=2008) & (treatgroup2==1|treat==0) & roadv_group==2  & lnNDVI_mn!=., a(county_num) rob cluster(county_num) 
	
	** Col.(3) High 
   xi: areg lnNDVI_mn ///
    P_x_T ///
	rain_m5 rain_m6 rain_m7 rain_m8 rain_m9 rain_m10 Tavg_m5 Tavg_m6 Tavg_m7 Tavg_m8 Tavg_m9 Tavg_m10 PSDI_m5 PSDI_m6 PSDI_m7 PSDI_m8 PSDI_m9 PSDI_m10 fiscal_ave ///
	i.year   ///
	if (year<=2013 & year>=2008) & (treatgroup2==1|treat==0) & roadv_group==3  & lnNDVI_mn!=., a(county_num) rob cluster(county_num) 

***Panel B. By NDVI in 2008： 
   ** Col.(1) low 
   xi: areg lnNDVI_mn ///
    P_x_T ///
	rain_m5 rain_m6 rain_m7 rain_m8 rain_m9 rain_m10 Tavg_m5 Tavg_m6 Tavg_m7 Tavg_m8 Tavg_m9 Tavg_m10 PSDI_m5 PSDI_m6 PSDI_m7 PSDI_m8 PSDI_m9 PSDI_m10 fiscal_ave ///
	i.year   ///
	if (year<=2013 & year>=2008) & (treatgroup2==1|treat==0) & ndvi08_group==1  & lnNDVI_mn!=., a(county_num) rob cluster(county_num) 
 
   ** Col.(2) medium  
   xi: areg lnNDVI_mn ///
    P_x_T ///
	rain_m5 rain_m6 rain_m7 rain_m8 rain_m9 rain_m10 Tavg_m5 Tavg_m6 Tavg_m7 Tavg_m8 Tavg_m9 Tavg_m10 PSDI_m5 PSDI_m6 PSDI_m7 PSDI_m8 PSDI_m9 PSDI_m10 fiscal_ave ///
	i.year   ///
	if (year<=2013 & year>=2008) & (treatgroup2==1|treat==0) & ndvi08_group==2  & lnNDVI_mn!=., a(county_num) rob cluster(county_num) 

   ** Col.(3) high  
   xi: areg lnNDVI_mn ///
    P_x_T ///
	rain_m5 rain_m6 rain_m7 rain_m8 rain_m9 rain_m10 Tavg_m5 Tavg_m6 Tavg_m7 Tavg_m8 Tavg_m9 Tavg_m10 PSDI_m5 PSDI_m6 PSDI_m7 PSDI_m8 PSDI_m9 PSDI_m10 fiscal_ave ///
	i.year   ///
	if (year<=2013 & year>=2008) & (treatgroup2==1|treat==0) & ndvi08_group==3  & lnNDVI_mn!=., a(county_num) rob cluster(county_num) 
	
***Panel C. By grassland type：
   ** Grassland type=meadow
	   xi: areg lnNDVI_mn ///
			 P_x_T ///
			 rain_m5 rain_m6 rain_m7 rain_m8 rain_m9 rain_m10 Tavg_m5 Tavg_m6 Tavg_m7 Tavg_m8 Tavg_m9 Tavg_m10 PSDI_m5 PSDI_m6 PSDI_m7 PSDI_m8 PSDI_m9 PSDI_m10  fiscal_ave ///
			 i.year ///
			 if (year<=2013 & year>=2008) & ((treatgroup2==1& grstype==1)|treat==0) ,a(county_num) rob cluster(county_num) 
			 
   ** Grassland type=grassland 
        xi: areg lnNDVI_mn ///
			 P_x_T ///
			 rain_m5 rain_m6 rain_m7 rain_m8 rain_m9 rain_m10 Tavg_m5 Tavg_m6 Tavg_m7 Tavg_m8 Tavg_m9 Tavg_m10 PSDI_m5 PSDI_m6 PSDI_m7 PSDI_m8 PSDI_m9 PSDI_m10  fiscal_ave ///
			 i.year ///
			 if (year<=2013 & year>=2008) & ((treatgroup2==1& grstype==2)|treat==0) ,a(county_num) rob cluster(county_num) 
			 
   ** Grassland type=desert 
		xi: areg lnNDVI_mn ///
			 P_x_T ///
			 rain_m5 rain_m6 rain_m7 rain_m8 rain_m9 rain_m10 Tavg_m5 Tavg_m6 Tavg_m7 Tavg_m8 Tavg_m9 Tavg_m10 PSDI_m5 PSDI_m6 PSDI_m7 PSDI_m8 PSDI_m9 PSDI_m10  fiscal_ave ///
			 i.year ///
			 if (year<=2013 & year>=2008) & ((treatgroup2==1& grstype==3)|treat==0) ,a(county_num) rob cluster(county_num) 
			 
   ** Grassland type=shub  
        xi: areg lnNDVI_mn ///
			 P_x_T ///
			 rain_m5 rain_m6 rain_m7 rain_m8 rain_m9 rain_m10 Tavg_m5 Tavg_m6 Tavg_m7 Tavg_m8 Tavg_m9 Tavg_m10 PSDI_m5 PSDI_m6 PSDI_m7 PSDI_m8 PSDI_m9 PSDI_m10  fiscal_ave ///
			 i.year ///
			 if (year<=2013 & year>=2008) & ((treatgroup2==1& grstype==4)|treat==0) ,a(county_num) rob cluster(county_num) 
			 
   ** Grassland type=Herbosa  
		xi: areg lnNDVI_mn ///
			 P_x_T ///
			 rain_m5 rain_m6 rain_m7 rain_m8 rain_m9 rain_m10 Tavg_m5 Tavg_m6 Tavg_m7 Tavg_m8 Tavg_m9 Tavg_m10 PSDI_m5 PSDI_m6 PSDI_m7 PSDI_m8 PSDI_m9 PSDI_m10  fiscal_ave ///
			 i.year ///
			 if (year<=2013 & year>=2008) & ((treatgroup2==1& grstype==5)|treat==0) ,a(county_num) rob cluster(county_num) 


			 
***** Supplementary Fiugre 1. Event study 
gen province_id_year=province_id*year
xi: areg lnNDVI_mn ///
     before5_P before4_P before3_P before2_P before1_P  after2_P after3_P after4_P after5_P ///
	 rain_m5 rain_m6 rain_m7 rain_m8 rain_m9 rain_m10 Tavg_m5 Tavg_m6 Tavg_m7 Tavg_m8 Tavg_m9 Tavg_m10 PSDI_m5 PSDI_m6 PSDI_m7 PSDI_m8 PSDI_m9 PSDI_m10  fiscal_ave ///
	 i.year ///
	 if (year<=2015 & year>=2006) & (treatgroup2==1|treat==0),a(county_num) rob cluster(province_id_year) 

	 
	 
***** Supplementary Table 1
// Note: these are obtained using data from 2006-2015,
// which is longer than the time period for the main analysis (2008-2013)

* column (1)
// row (1)
sum lnNDVI_mn if (P==1) & (year>=2006 & year<=2010)
// row (2)
sum lnNDVI_mn if (P==1) & (year<=2015 & year>=2011)
// row (3)
reg lnNDVI_mn T if (year>=2006 & year<=2015) & P==1,rob cluster(county_num)

* column (2)
* row (1)
sum lnNDVI_mn if P==0 & (year>=2006 & year<=2010)
* row (2)
sum lnNDVI_mn if P==0 & (year<=2015 & year>=2011)
* row (3)
reg lnNDVI_mn T if (year>=2006 & year<=2015) & P==0,rob cluster(county_num)
	
* colume (3)
* row (1)
reg lnNDVI_mn P if (year>=2006 & year<=2010) & (P==1|P==0),rob cluster(county_num)
* row (2)
reg lnNDVI_mn P if (year>=2011 & year<=2015) & (P==1|P==0),rob cluster(county_num)
* row (3)
xi: reg lnNDVI_mn P T P_x_T if (year<=2015 & year>=2006),rob cluster(county_num)




*****Supplementary Table 2. Robustness Checks: different treatments and controls.
		*** Col.(1): Inner Mongolia 
		
		xi: areg lnNDVI_mn ///
			 P_x_T ///
			 rain_m5 rain_m6 rain_m7 rain_m8 rain_m9 rain_m10 Tavg_m5 Tavg_m6 Tavg_m7 Tavg_m8 Tavg_m9 Tavg_m10 PSDI_m5 PSDI_m6 PSDI_m7 PSDI_m8 PSDI_m9 PSDI_m10  fiscal_ave ///
			 i.year ///
			 if (year<=2013 & year>=2008) & (treatgroup3==1|treat==0),a(county_num) rob cluster(county_num) 
			 
		*** Col.(2): Bordered 
		xi: areg lnNDVI_mn ///
			 P_x_T ///
			 rain_m5 rain_m6 rain_m7 rain_m8 rain_m9 rain_m10 Tavg_m5 Tavg_m6 Tavg_m7 Tavg_m8 Tavg_m9 Tavg_m10 PSDI_m5 PSDI_m6 PSDI_m7 PSDI_m8 PSDI_m9 PSDI_m10  fiscal_ave ///
			 i.year ///
			 if (year<=2013 & year>=2008) & (treatgroup3==1|treat==0)& border==1,a(county_num) rob cluster(county_num) 
			 

***** Supplementary Table 3. Robustness check
xi: areg lnNDVI_mn ///
     P_x_T ///
	 rain_m5 rain_m6 rain_m7 rain_m8 rain_m9 rain_m10 Tavg_m5 Tavg_m6 Tavg_m7 Tavg_m8 Tavg_m9 Tavg_m10 PSDI_m5 PSDI_m6 PSDI_m7 PSDI_m8 PSDI_m9 PSDI_m10  fiscal_ave ///
	 i.year ///
	 if (year<=2011 & year>=2010) & (treatgroup2==1|treat==0),a(county_num) rob cluster(county_num) 

 xi: areg lnNDVI_mn ///
     P_x_T ///
	 rain_m5 rain_m6 rain_m7 rain_m8 rain_m9 rain_m10 Tavg_m5 Tavg_m6 Tavg_m7 Tavg_m8 Tavg_m9 Tavg_m10 PSDI_m5 PSDI_m6 PSDI_m7 PSDI_m8 PSDI_m9 PSDI_m10  fiscal_ave ///
	 i.year ///
	 if (year<=2012 & year>=2010) & (treatgroup2==1|treat==0),a(county_num) rob cluster(county_num) 

  xi: areg lnNDVI_mn ///
     P_x_T ///
	 rain_m5 rain_m6 rain_m7 rain_m8 rain_m9 rain_m10 Tavg_m5 Tavg_m6 Tavg_m7 Tavg_m8 Tavg_m9 Tavg_m10 PSDI_m5 PSDI_m6 PSDI_m7 PSDI_m8 PSDI_m9 PSDI_m10  fiscal_ave ///
	 i.year ///
	 if (year<=2013 & year>=2010) & (treatgroup2==1|treat==0),a(county_num) rob cluster(county_num) 

  xi: areg lnNDVI_mn ///
     P_x_T ///
	 rain_m5 rain_m6 rain_m7 rain_m8 rain_m9 rain_m10 Tavg_m5 Tavg_m6 Tavg_m7 Tavg_m8 Tavg_m9 Tavg_m10 PSDI_m5 PSDI_m6 PSDI_m7 PSDI_m8 PSDI_m9 PSDI_m10  fiscal_ave ///
	 i.year ///
	 if (year<=2014 & year>=2010) & (treatgroup2==1|treat==0),a(county_num) rob cluster(county_num) 

   xi: areg lnNDVI_mn ///
     P_x_T ///
	 rain_m5 rain_m6 rain_m7 rain_m8 rain_m9 rain_m10 Tavg_m5 Tavg_m6 Tavg_m7 Tavg_m8 Tavg_m9 Tavg_m10 PSDI_m5 PSDI_m6 PSDI_m7 PSDI_m8 PSDI_m9 PSDI_m10  fiscal_ave ///
	 i.year ///
	 if (year<=2015 & year>=2010) & (treatgroup2==1|treat==0),a(county_num) rob cluster(county_num) 
	 

***** Supplementary Table 5-7. Summary statistics 
		gen treat_11=0 if treat==0 
		replace treat_11=1 if treatgroup2==1 

		gen treat_22=0 if treat==0 
		replace treat_22=1 if treatgroup3==1 
		
		gen treat_33=0 if treat==0 & border==1
		replace treat_33=1 if treatgroup3==1 & border==1
		
				
		sum rain_m5-rain_m10 Tavg_m5-Tavg_m10 PSDI_m5-PSDI_m10 fiscal_ave if year==2010 & treat_11==0 &lnNDVI_mn!=.
		sum rain_m5-rain_m10 Tavg_m5-Tavg_m10 PSDI_m5-PSDI_m10 fiscal_ave if year==2010 & treat_11==1 &lnNDVI_mn!=.
		
		
		sum rain_m5-rain_m10 Tavg_m5-Tavg_m10 PSDI_m5-PSDI_m10 fiscal_ave if year==2010 & treat_22==0 &lnNDVI_mn!=.
		sum rain_m5-rain_m10 Tavg_m5-Tavg_m10 PSDI_m5-PSDI_m10 fiscal_ave if year==2010 & treat_22==1 &lnNDVI_mn!=.
		
		sum rain_m5-rain_m10 Tavg_m5-Tavg_m10 PSDI_m5-PSDI_m10 fiscal_ave if year==2010 & treat_33==0 &lnNDVI_mn!=.
		sum rain_m5-rain_m10 Tavg_m5-Tavg_m10 PSDI_m5-PSDI_m10 fiscal_ave if year==2010 & treat_33==1 &lnNDVI_mn!=.


/* FOR AUTHOR THEMSELVES
keep lnNDVI_mn P_x_T rain_m5 rain_m6 rain_m7 rain_m8 rain_m9 rain_m10 Tavg_m5 Tavg_m6 Tavg_m7 Tavg_m8 Tavg_m9 Tavg_m10 PSDI_m5 PSDI_m6 PSDI_m7 PSDI_m8 PSDI_m9 PSDI_m10  fiscal_ave /// 
year province_id county_num  treatgroup2 treatgroup3 border treat lncattle lnsheep roadv_group ndvi08_group  grstype ///
before5_P before4_P before3_P before2_P before1_P  after2_P after3_P after4_P after5_P 


label var province_id "Province ID"
label var year "Year"
label var county_num "county ID"

label var PSDI_m5 "PSDI drought index for May"
label var PSDI_m6 "PSDI drought index for June"
label var PSDI_m7 "PSDI drought index for July"
label var PSDI_m8 "PSDI drought index for August"
label var PSDI_m9 "PSDI drought index for September"
label var PSDI_m10 "PSDI drought index for October"

label var Tavg_m5 "mean temperature for May"
label var Tavg_m6 "mean temperature for June"
label var Tavg_m7 "mean temperature for July"
label var Tavg_m8 "mean temperature for August"
label var Tavg_m9 "mean temperature for September"
label var Tavg_m10 "mean temperature for October"

label var rain_m5 "cumulative precipitation for May"
label var rain_m6 "cumulative precipitation for June"
label var rain_m7 "cumulative precipitation for July"
label var rain_m8 "cumulative precipitation for August"
label var rain_m9 "cumulative precipitation for September"
label var rain_m10 "cumulative precipitation for October"

label var fiscal_ave "Per-capita county fiscal income (yuan)"

label var treat "treat=0 if control group; treat=1 if treatment group"

label var border "=1 if bordered counties"

label var P_x_T "DID key variable; =1 if treated after treatment period; otherwise=0"

label var treatgroup2 "use the second group as the treatment group; =1 if treatment group,=0 otherwise"
/* The treatment group (P = 1) includes the counties in five North and Northwestern program provinces that were covered in GECP-I, i.e., Xinjiang, Qinghai, Gansu, Ningxia, and Inner Mongolia. The control group (P = 0) includes the counties in five North and Northeastern provinces that were not covered by GECP-I until 2016, i.e., Shanxi, Hebei, Liaoning, Jilin, and Heilongjiang.*/

label var treatgroup3 "treatment grous are bordered treatment provinces"

label var lnNDVI_mn "log of NDVI"

label var lncattle "log of cattle inventory"
label var lnsheep  "log of sheep inventory"

label var grstype "grassland type 1=meadow; 2=grassland; 3=desert; 4=shrub; 5=herbosa"

label var roadv_group "road intensity groups: 1=low; 2=medium; 3=high"

label var ndvi08_group "NDVI in 2008 by goups: 1=low; 2=medium; 3=high"

label var before5_P "=1 if 5 years before the program started & treated" 
label var before4_P "=1 if 4 years before the program started & treated" 
label var before3_P "=1 if 3 years before the program started & treated" 
label var before2_P "=1 if 2 years before the program started & treated" 
label var before1_P "=1 if 1 years before the program started & treated" 
label var after2_P "=1 if the second year after the program started & treated" 
label var after3_P "=1 if the third year after the program started & treated" 
label var after4_P "=1 if the fourth year after the program started & treated" 
label var after5_P"=1 if the fifth year after the program started & treated" 
*/

